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ABSTRACT 

We use a fully 3-dimensional thermal evolution model to examine the effects 
of a non-uniform surface albedo on the subsurface thermal structure of comets. 
Surface albedo markings cast "thermal shadows", with strong lateral thermal 
gradients. Corresponding compositional gradients can be strong, especially if the 
crystallization of amorphous water ice is triggered in the hottest regions. We show 
that the spatial extent of the structure depends mainly on the obliquity, ther- 
mal conductivity and heliocentric distance. In some circumstances, subsurface 
structure caused by the thermal shadows of surface features can be maintained 
for more than 10 Myr, the median transport time from the Kuiper Belt to the 
inner solar system. Non-uniform compositional structure can be an evolutionary 
product and does not necessarily imply that comets consist of building blocks 
accumulated in different regions of the protoplanetary disk. 

Subject headings: Comets: general, methods: numerical 
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Introduction 



Solar system comets are currently stored in two main reservoirs, namely the Oort 
Cloud and the Kuiper Belt, having different dynamical histories and physical properties. 
The bodies contained in these reservoirs can be scattered to the inner solar system by 
various gravitational processes. The Scattered Disk component of the Kuiper Belt is widely 
believed to be the source of Jupiter Family Comets (hereafter JFCs, Volk fc Malhotra]|2008 ). 
In particular, a gravitational cascade might exist between the three distinct populations: 



scattered disk objects. Centaurs and JFCs (Levison & Duncan 1997 Tiscareno & Malhotra 



2003). 



Comets are believed to contain some of the best-preserved material from the formation 
of our planetary system. Cometary composition should reflect the location at which the 
material -ices and dust- formed in the protoplanetary disk. A gradient in composition 
might reasonably be expected, distinguishing comets formed at high temperatures close to 
the proto-Sun, from those formed in an environment dominated by interstellar chemistry at 
large heliocentric distances. This simple picture is complicated by potential radial mixing 
inside the protoplanetary disk. On a micro-scale, the co-existence of crystalline silicates 
(formed at temperatures > 10^ K) with cometary ice (accreted at temperatures < 50 K) 
provides direct evidence for radial mixing dCampins fc Ryan||l989| |Ishii et al.||2008| . On 
a macro-scale, cometesimals formed at different heliocentric distances and temperatures 
might have been scattered and later aggregated to form individual cometary nuclei. For 
this reason, one key question in cometary science is whether the comet nuclei observed 
today are internally homogeneous or heterogeneous in composition. Observationally, both 
types of chemical structure have been reported among the JFCs. For example outgassing 



from comets 9P/Tempel 1 (Mumma et al. 2005 Feaga et al. 2007) and 103P/Hartley 
2 (A'Hearn et al. 2011) appears to be compositionally non-uniform, whereas comet 
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73P/Schwassmann-Wachmann 3 is found to be uniform (DelloRusso et al. 2007 Kobayashi 



et al. 2007) 



However, several post-accretion processes could alter the primordial compositions of 
comets. Their surface layers might be chemically stratified by solar wind and cosmic ray 
irradiation during their 4.5 Gyr residence in the Kuiper Belt and the Oort Cloud reservoirs 



(Cooper et al. 2003). In addition, their internal structures and compositions could be 



locally modified by heat absorbed at the surface from the Sun (for example, images of the 



nucleus of comet 103P/Hartley 2 show local albedo variations up to a factor of 4 (A'Hearn 



et al. 2011). The diversity in composition observed in comets (A'Hearn et al. 1995) could 



consequently be the result of a combination of these three effects: different formation 
environments, chemical evolution during multi-billion year storage in the source reservoirs, 
and recent thermal processing by absorbed sunlight once injected into the planetary region. 

In this paper, we use a new thermal evolution model that has been developed to allow 



fast, accurate computation of the 3-dimensional heat transport problem ( Guilbert-Lepoutre 



et al. 2011). The speed of the model allows us to accurately calculate the effects of 



non-uniform surface albedo on the internal structure of a cometary nucleus and to explore 
the influence of orbital and thermophysical parameters. The model and assumptions are 
presented in Section 2. The results for different configurations are presented in Section 3 
and discussed in Section 4. 



2. Thermal evolution model 

2.1. Main equations 

Ours is a "toy-model" in which the parameters of the nucleus are idealized in order to 
make the problem tractable, and no attempt is made to model any particular real nucleus. 
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The body is assumed to be initially a sphere made of a porous mixture of ice and dust 



uniformly distributed within the icy matrix. Jewitt (2009) and Meech et al. (2009) report 



observational, albeit indirect, evidence consistent with the presence of amorphous water ice 
in comets. We therefore assume the ice is initially amorphous. The model we use evaluates 
the temperature distribution inside the body by taking into account three dimensional 
heat fluxes, and includes energy release from the crystallization of amorphous ice. This 



3-dimensional model is fully described in Guilbert-Lepoutre et al. (2011). We here give only 



an outline. The heat conduction equation to be solved is the following: 

dT 



pbuikC -1— + V(-K T) 
at 



crysti 



where T [K] is the temperature distribution to be determined, pi^ik [kg m~^] the object's 
bulk density, c [J kg^^ K^^] the material heat capacity, n [W m~^ K~^] its effective thermal 
conductivity (parameters described next section), and Qcryst [W m^^] the internal power 
production per unit volume due to the amorphous-crystalline phase transition. The latter 
is described by: 

Qcryst = A(T) Pa Hac, (2) 

with Pa [kg m~^] the amorphous water ice bulk density. The phase transition releases a 
latent heat Hac = 9x10"^ J kg^^ ( Klinger|[l981 ), at a rate measured by Schmitt et al. (1989): 



A(T) = 1.05 X 10^3 e-5370/T 



(3) 



Boundary conditions are considered both at the surface and in the center of the object. 
Several thermal processes are considered to evaluate the thermal balance for each point on 
the surface: 



solar illumination described by (1 — cos^, with A the Bond albedo, Sq the solar 
constant, dn the object's heliocentric distance, and C, < 90° the local zenith angle. 
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- thermal emission eaT^, with e material emissivity, a the Stefan-Boltzmann constant 
and T the surface temperature. 

- lateral and radial heat fluxes. 

The heat diffusion equation ([T]) is expanded in spherical coordinates: 

PbulkC dT (2 ^ ldK\dT d^T 1 ^ ^ _ Q,ryst 



K dt \r ' K dr J dr dr^ ^2 k ' 

with Ag the angular Laplacian operator. As spherical harmonics allow a simple and 
natural expression of the temperature over a regular spherical grid, we introduce them to 

describe the temperature distribution: 

00 I 

T = Y.Y. T'-it,r)yUO,^). (5) 

1=0 m=-l 

This sum is exact as long as the degree I goes to inflnity. This cannot be reached in practise 

and the sum is thus cut to a maximum degree Imax- We then introduce the expansion in the 

heat diffusion equation: 

PC aT'- f2 , ldK\ dT^"^ m + 1) ,^ 

+ ——o — T = Q (6) 



K dt \r Kdr J dr dr^ r^ 
with Q^'^yst = Qcryst Sifl 5mfl, S being the Kronecker function. We therefore obtain 
{imax + 1)^ equations of T^"'{t,r), instead of one single 3-dimensional equation for T. These 
ID equations are solved using a Crank-Nicholson numerical scheme, which is a stable 
implicit technique. 

The boundary condition at the surface is given by the thermal balance evaluated for 
each point of the surface, expanded into the basis of spherical harmonics: 

Jo Jo 



The sampling theorem developed by DriscoU & Healy (1994) is used to derive these 



coefficients. Denoting by the number of points in one direction of the equally sampled 
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surface grid, the boundary conditions Tjj^j are computed as: 

N-l N-l 

TiZf = E E «f ' ^^-rf{0j, V>k) yimiOj, V>k), (8) 
j=0 k=0 

with 6j = (fk = ^ the grid point coordinates, aj^"^ a coefficient that accounts for the 
over-samphng near the poles, and Tg^rf [K] the equihbrium temperature at each point of 
the surface. The number of points in one angular direction is chosen so as to minimize 
the discretization errors and the computational load. In the center, the boundary condition 
is simpler and depends only on r: 

dr dr 



2.2. Thermophysical properties 

The bulk density is related to the porosity of the solid matrix ip by: 

= (1 - ^) f ^ + ^ V\ (10) 

V PH20 pd J 

with Xu^o and the mass fractions of water ice and dust respectively, PH20 and pd [kg 
m~'^] the densities of water ice and dust respectively. We assume that the object is made of 
a material with a dust to water ice mass ratio XdjXu^o = 1, a porosity i\} = 30%, and a 
bulk density pbuik = 1 g cm"^. The heat capacity of the mixture is obtained by computing 
the average of the values weighted by the mass fraction of each component 

C = XH2OCH2O + XdCd, (11) 

with and the mass fraction of water ice and dust, and CH2O and c^ [J kg"^ K^^] the 

heat capacities of each component. The numerical values used in this work can be found in 
Table [H 
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We evaluate the thermal conductivity by considering the material as made of two 
phases, the empty pores with a thermal conductivity and the solid matrix with a thermal 
conductivity Hg. Within the empty pores the heat is transferred through thermal radiation 
for which the effective conductivity is 



Kp = ArpSaT 



(12) 



with Tp =1 /zm the average pore radius, e=0.9 the medium emissivity, a the Stefan- 



Boltzmann constant and T [K] the temperature (Huebner et al. 2006). The solid 
matrix thermal conductivity, k^, is computed as the average of each component thermal 
conductivity (see Table [T]), weighted by its volume fraction: 



XH2O [(1 - Xcr)lia + XctI^ct] + Xdl^d, 



(13) 



with XH2O and Xd the volume fractions of water ice and dust respectively, and X^r the mass 
fraction of crystalline water ice. We also consider a Hertz factor h with a fixed value of 0.1 



to account for the granular structure of the solid (Huebner et al. 2006). We finally use the 



Russel formula (Russel 1935) to calculate a correction factor which should be applied to 



Ks, to account for the effects of porosity (Espinasse et al. 1991 Coradini et al. 1997 Orosei 



et al. 



1999). It depends on the porosity and the ratio / 

^2/3j ^ (1 _ ^2/3) 



as 



^ ^ - ^2/3 + 1 _ ^2/3(^1/3 _ 

The material effective thermal conductivity is consequently: 



(14) 



h 



(15) 



2.3. The object, the orbit and the albedo patch at the surface 

As mentioned in Section 2.1, the object is assumed to be a sphere, made of a porous 
matrix of amorphous water ice and dust. The thermal properties of such a mixture have 



- 9 - 



been presented in the previous section. We also assume a radius R oi 2 km for the object, 



which is a typical comet radius (A'Hearn et al. 1995). Following the idea that there 



might be a dynamical cascade from the Kuiper Belt to the Centaurs to the JFCs, due to 
gravitational interactions with giant planets, we assume that the body enters the inner 
solar system on a Centaur-like orbit. The rotation period of the object is considered to be 



10 hrs, typical of Centaurs and KBOs (Sheppard et al. 2008). 



Finally, we assume that the initial object has a non-uniform Bond albedo at the 
surface. We consider a surface of 10% Bond albedo (Agurf), with a patch of 60% Bond 
albedo (Apatch), corresponding to a fresh ice/frost region on an otherwise dirty ice or 
refractory surface. While most TNOs and Centaurs have dark surfaces, some have high 



average albedos due to the presence of surface ice (Stansberry et al. 2008). The patch is 



initially positioned between latitudes +22.5° and -22.5°, and between longitudes 0° and 
180°. This corresponds to l/8th of the overall surface, leading to a spherical average albedo 
of about 16%. The parameters of the object, the orbit and the patch are summarized in 
Table |2} For each simulation we also consider a reference case in which all the parameters 
are the same except that there is no albedo spot at the surface. 



3. Results 

In the interests of brevity, we focus the discussion on six cases as summarized in Table 
[3] (c.f. Figures [1] to [7]) . Cases A, D, E and F illustrate the influence of the heliocentric 
distance through variations of the semimajor axis and the eccentricity. Case B illustrates 
the effects of the obliquity, and Case C illustrates the influence of the material thermal 
conductivity. We performed additional simulations to explore the effects of the size and 
position of the patch or the albedo difference between the patch and the surface. 



3.1. Effects of non-uniform albedo 



The spatially varying albedo induces a diurnally and annually modulated heat wave 
in the nucleus, with a peculiar shape that produces lateral subsurface thermal gradients. 
These effects are illustrated by Figsjl] and [2] for Case A. Figure [T] shows the evolution of 
the temperature beneath a given point on the equator with or without (reference case) the 
albedo spot. The temperature difference at the surface caused by the spot is about 30 K. 
Figure [2] shows that the region beneath the albedo patch remains 20 to 30 K cooler than 
in the reference case, creating a thermal shadow that appears very quickly when the object 
enters the inner solar system. On both figures, the black line labeled H20cr stands as a 
limit beyond which amorphous water ice has been crystallized (in the hottest regions). 

We expect that compositional gradients would also develop, following these lateral 
thermal gradients, because water ice crystallization and volatile sublimation are strongly 
temperature dependent. Specifically, cold regions in the thermal shadows of surface albedo 
features should be enhanced in volatiles relative to neighboring unshadowed regions, for 
two reasons. First, uncrystallized ice in a local cold spot will retain its full complement 
of trapped volatiles, while these volatile species will have been liberated from surrounding 
crystallized ice. Second, the thermal shadows may act as cold-traps in which volatile 
abundances are enhanced even further by the diffusive migration of molecules from 
crystallized ice adjacent to uncrystallized ice. 

We consider these points for Case A, in which crystallization is triggered by insolation 
both in the reference and non-uniform-albedo models, while the region located under the 
albedo spot remains cool enough to prevent crystallization (Fig|2]). Volatiles trapped in the 



amorphous matrix will be released upon crystallization (Bar- Nun et al. 1985 Laufer et al. 



1987 Notesco & Bar-Nun 1996 Bar-Nun & Owen 1998 Notesco et al. 2003) and travel in 



a free molecular flow (Knudsen flow, which is typical for comets, see Huebner et al. (2006)). 
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The diffusion coefficient of gases released upon crystallization is given by Prialnik (1992): 

4 



D 



3 ^ 



(16) 



with V 



being the mean thermal gas velocity {ks is the Boltzmann constant, 



T [K] the temperature and m [g] the molecule mass), and Kp a length coefficient which 
characterizes the porous material. For a medium made of randomly packed spheres with a 
resulting porosity ip and a pore radius rp, Kp = ^^J^y/s (Prialnik 



1992). The pore radius 



is very uncertain. We here use rp=10~^ m but values an order of magnitude smaller are 
possible. Substituting ip=30% gives Kp = 1.85 x 10^^ m. The diffusion coefficients of 
CO and CO2 are Dco=7.11x 10"^ m^ s"^ and i5c'02=5-67x 10"^ m^ s~^ respectively, with 
T=110 K the maximum temperature reached within the crystalline regions. 



The gas diffusion length is given hj i = I^Dt with t [s] the time. For example, 
in a typical nucleus rotation period of ~10 hr, the molecules can flow through i ~3 m 
of the porous material given the above diffusion coefficients. In one orbit (18.5 yr for 
a=7 AU), CO and CO2 molecules can flow through i ~400 m. Each thermal shadow is thus 
surrounded by a layer of radial extent I from which liberated volatile molecules might be 
trapped. The trapping will not be perfectly efficient, because liberated molecules can also 
migrate to the free surface of the nucleus and escape, or move down the thermal gradient 
into the cometary interior and re-freeze. Nevertheless, we expect that migration into 
thermal shadow cold-traps will produce a preferred spatial scale for volatile segregation, 
with the volatile enhancement being strongest for albedo spots having size comparable to 
the diffusion length. For larger albedo spots, we expect that migrating volatiles will be 
trapped in a rim having thickness comparable to i. 

On average, a Centaur can spend ~10 Myr on its orbit before either leaving the solar 



system or becoming a JFC, due to gravitational interactions with giant planets (Tiscareno & 



Malhotra 2003 Horner et al. 2004). We performed simulations over 10 Myr, which showed 
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that the subsurface cold plug can persist. After 10 Myr on the Centaur orbit, the thermal 
shadow produced by the presence of the higher albedo spot at the surface reaches about 
400 m deep provided the nucleus spin vector remains constant over this period. Lateral 
heat fluxes only begin to erase the lateral thermal gradients, which are still important, in 
particular close to the surface where the temperature difference is about 20 to 30 K. 



3.2. Effect of obliquity 

We found that obliquity is the parameter most affecting the formation of a thermal 
shadow. At non-zero obliquity, the variations in the subsolar point latitude introduce 
asymmetry in the propagation of the heat wave, as illustrated by Case B in Figjs] (0 = 
20°). The temperature distribution with a non-zero obliquity varies across the orbit due 
to the variations of the subsolar point latitude, in addition to variations attributable to 
the propagation of the heat wave (which have been illustrated previously). Figure [s] thus 
corresponds to a snapshot of the distribution, in which the subsolar point is located at the 
equator, and is moving southward. 

After 10 Myr spent in the inner solar system, the temperature distribution for non-zero 
obliquity tends to become more uniform as the effects of the latitudinal movements of 
the subsolar point are averaged out. While the temperature difference between the poles 
and the equator is almost 40 K for O = 0°, it is only 5 K with O = 45° after 10 Myr. 
Consequently, for high obliquities (typically larger than 45°), the thermal shadow disappears 
with time, even if it appeared during the first few orbits. In the extreme case of G = 
90°, no thermal shadow ever appears despite the presence of temporary strong lateral 
thermal gradients. The radial propagation of the heat wave and lateral heat fluxes erase 
the potential compositional gradients already during the first orbit. Still we find that a 
subsurface thermal shadow can be maintained over 10 Myr on a Centaur-like orbit if the 
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obliquity is low (< 30° for substantial volatile enhancement effect). 

3.3. Variations of other parameters 

The thermal conductivity affects mainly the radial extent of the cold plug, since it 
controls the efficiency of the heat transfer in the material, as illustrated by Case C and 
Figjlj In this case, we considered a thermal conductivity ten times lower than in the other 
cases. Interestingly, a low thermal conductivity implies that lateral heat fluxes are very 
ineffective in erasing any plug thus produced, which can survive for more than 10 Myr even 
if their lateral extent is initially small. 

The heliocentric distance is also an important factor. With increasing heliocentric 
distance, the steepness of the lateral thermal gradients decreases. In addition, as shown by 
FigjSjfor Case D, crystallization might not be reached in the hottest regions, thus limiting 
the volatile enhancement in the thermal shadows of surface features. On the contrary, 
a smaller heliocentric distance will induce higher surface temperatures. In Case E, the 
crystallization threshold is reached also under the albedo spot (Fig|6]). Over the time of 
residence in the giant planets region, the non-uniform structure could still exist. The most 
volatile enriched region would nonetheless be located a few meters to a few tens of meters 
deep. We found that the orbital eccentricity has very limited effect (Case F). Its influence 
is restricted to the amount of energy to be transferred to the subsurface (see Figj?]), which 
varies around the orbit. This merely impacts the shape or extent of the temperature 
distribution in the case we show in FigjT} The orbit that we considered has a semi-major 
axis of 15 AU and an eccentricity of 0.2, resulting in a perihelion distance of 12 AU, the 
same distance as in Case D for a circular orbit with a=12 AU. Larger eccentricities could 
nonetheless have more influence, as the energy provided close to perihelion could potentially 
trigger the crystallization of amorphous water ice. 
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The surface of a comet is impacted by a large variety of energetic particles, which 



might produce an irradiation crust (StrazzuUa et al. 1991 Hudson et al. 2008). The first 
few centimeters are the most affected, but high energy particles might penetrate up to a 
meter beneath the surface. The thermal properties of such a crust are not constrained yet, 
but recent laboratory experiments on porous dust aggregates indicate thermal conductivity 



between 10~ and 10~ Wm~ K~ (Krause et al. 2011). This is very similar to the thermal 
conductivity of the material we are considering. Therefore, this crust would not prevent 
the progression of the heat wave toward the center of the object. It would have a damping 
effect, moving the interesting boundary a meter deeper, and inducing lower temperatures. 
The extent of non-uniform structures could be limited in this case, but would still exist, 
especially for objects orbiting close to the Sun (as in Case E). 



4. Discussion 

Our simulations show that non-uniform surface albedo creates thermal gradients in the 
subsurface layers, which can produce a long-lived non-uniform subsurface structure. The 
emergence of compositional non-uniformity depends on the albedo difference between the 
surface and the patch, rather than the albedo itself. Local albedo variations of a factor of 



four exist on comet 103P/Hartley 2 (A'Hearn et al. 2011), suggesting that thermal shadow 



effects could be very strong. Although we considered a large albedo spot, smaller scale 
surface features would have the same impact on the subsurface, since with low thermal 
conductivities, lateral heat transfers are quite inefficient in erasing the cold plugs. The 
overall albedo could remain very low, as the variation attributable to the spot would 
be limited. Real comet nuclei have complex shapes and surface features, as revealed by 
spacecraft observations like for comets 9P/Tempel 1 or 103P/Hartley 2 (Deep Impact and 
EPOXI missions respectively). The effects of surface topographic features like craters would 
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mimic those of albedo features, since ultimately the important factor is the heat transferred 
to the subsurface. Consequently, the presence of craters, topographic features, boulders 
or any other source of shadowing can generate subsurface lateral thermal gradients and 
non-uniform compositions. 

If the crystallization threshold is locally reached, the thermal shadow could become 
strongly enriched in volatiles, while the surrounding crystallized ice would be depleted in 
volatiles. In our cases, the super-volatile CO would most likely escape the body, while a 
less volatile compound such as CO2 could recondense in the thermal shadow. The diffusion 
length is independent of the patch size, and depends mainly on the material thermal and 
structural properties. Time is also an important factor. If the size of the spot is comparable 
to or smaller than the diffusion length, the volatile abundance enhancement of the cold plug 
could be extremely important. Such features have been reported to have a scale of tens to 



hundreds of meters (Mumma et al. 1993; Weissman et al. 2004 A'Hearn et al. 2011), which 



is very similar to the diffusion lengths considered here. Consequently, we can predict that 
the scale of the enriched regions would range from a few meters to a few hundred meters, 
depending on i) the size of the patch, as the size of the enriched region cannot be larger 
than the cold plug, ii) the material properties and iii) the time the molecules had to flow 
in the porous medium. Nonetheless, the gas phase is not accounted for nor modeled in our 
simulations. The effects of such a gas phase could strongly modify our results, as it can 
locally affect the thermal conductivity, the porosity or pore sizes. In addition, instabilities 
caused by pressure build-up could develop and blow up some surface layers. There are, 
however, too many unknowns to meaningfully model all these processes. 

The compositional variations produced in the subsurface can be sustained until the 
body becomes a JFC. In this case, the resulting internal composition would be strongly 
non-uniform, and the cometary activity would be generated through jets. These would 
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be produced as the object orbits closer to the Sun, where insolation can finally trigger 
the sublimation of volatiles which were concentrated in the thermal shadows, and/or the 
crystallization of these regions. Consequently, we expect that non-uniform thermal and 
compositional structure should be common. Identification of compositional differences in 
a single nucleus does not necessarily imply that the comets were built from cometesimals 
formed at different heliocentric distances, with distinct compositions. 



5. Summary 

Fast, 3-dimensional thermal evolution simulations show: 

1. Non-uniform surface albedos on comets can generate long-lived thermal shadows in 
the immediate sub-surface regions. 

2. Temperature-sensitive processes (including sublimation and crystallization) proceed 
at different rates inside and outside the thermal shadows, leading to the development 
of volatile-enhanced shadow cold-traps. 

3. Compositional gradients caused by thermal shadows should be most pronounced for 
albedo spot sizes comparable to the diffusion length (typically from a few meters to a 
few hundred meters for structures growing on the rotational and orbital timescales). 

4. Under some circumstances, subsurface temperature structure can be preserved for 
the mean hfetime of a Centaur (10 Myr) before the object becomes a Jupiter Family 
Comet. 

5. Observations of jets and non-uniform compositions in cometary nuclei do not 
necessarily imply an initially non-uniform composition. 
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Table 1: Heat capacities and thermal conductivities of different components in the material 
mixture. 



Param. 


Value 


Unit 




Ref. 


CH2O 


7.49T + 90 


Jkg-^K- 


1 


G&S36 


Cd 


1200 


Jkg-iR- 


1 


E&S83 




2.34 10-'^T + 21 


] 10-2 Wm-^K 


-1 


K180 


Kcr 


567 / T 


Wm-^K 


-1 


K180 




4.2 


Wm-^K 


-1 


E&S83 


Cinit 


760 


Jkg-^K- 


1 






6.17 10^2 


Wm-^K 


-1 





Note. — Heat capacities: cjj^o for water ice and Cd for dust respectively, thermal conductivities: Ka and 
Kcr for amorphous and crystalline water ice respectively, and for dust. Cinit and Kinit correspond to the 
initial values of the heat capacity and thermal conductivity in the simulations. 



References. — G&S36: Giauque k Stout (1936), E&S83: Ellsworth fc Schubert (1983), K180: Klinger 



(1980). 
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Table 2: Initial values for various parameters relative to the object, its orbit and the non- 
uniform surface albedo, considering Case A. 



Par am. 


R 


Pbulk 




'^patch 


a* 


a 


e e 






Unit 


km 


g cm""^ 








AU 





hrs 


yrs 


Value 


2 


1 


10% 


60% 


12.5% 


7 





10 


18.5 



Note. — 



area ratio between the spot and the overall surface. 
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Table 3: Parameters for the different cases illustrated, with the corresponding figures. 



Uase 


[VVni IV 


-11 
I 


a 

TATTl 
[AUJ 


e 


e Fi^ 

[1 


5. Description 


A 


6.17 10" 


_2 


7 








1 


i radial evolution with time 














2 


T distrib. after one orbital per. 


B 


6.17 10- 


-2 


7 





20 


3 


T distrib. after one orbital per. 


C 


6.17 10" 


3 


7 








4 


T distrib. after one orbital per. 


D 


6.17 10- 


^2 


12 








5 


T distrib. after one orbital per. 


E 


6.17 10- 


-2 


5 








6 


T distrib. after one orbital per. 


F 


6.17 10" 


^2 


15 


0.2 





7 


T radial evolution with time 
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Fig. 1. — Case A. Radial evolution of the temperature under a given point on the equator: 
reference model (top panel) and inside the high albedo spot (bottom panel). The water ice 
crystallization boundary is delineated by the black line. No ice crystallizes under the spot. 



-22- 



T 



cn 
CD 



D 



CD 
"D 

CD 




2 4 6 8 
Depth [m] 



-50 







2 4 6 8 

Depth [m] 







1 00 


■ 


90 




80 




70 




60 




50 


1 


40 





1 2 




12 



T 



1 00 


■ 


90 




80 




70 




60 




50 


1 


40 





Fig. 2. — Case A. Temperature distributions along a meridian: reference model (top panel) 
and with a surface albedo spot located between latitudes +22.5° and -22.5° (bottom panel), 
after one orbital period. The water ice crystallization boundary is delineated by the black 
line. 
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Fig. 3. — Case B. Same as Fig|2] showing asymmetries produced by obliquity 20°. The 
subsolar point is located at the equator moving southward. 
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Fig. 4. — Case C. Same as Fig|2]with a thermal conductivity 10 times lower than in Case 
A, showing the reduced spatial scale of the thermal structure. 
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Fig. 5. — Case D. Same as Fig|2]with a=12 AU. The crystallization threshold is not reached 
at this larger distance. 
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Case E. Same as Figj2]with a=5 AU. 
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Fig. 7. — Case F. Same as FiglTlfor an eccentric orbit (a=15 AU, e=0 
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